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Abstract 

This  paper  summarizes  the  current  effort  at  the  Air  Force  Research  Laboratory  at  Edwards  AFB  in 
developing  a  general  hybrid  model  for  the  studies  of  multi-scale,  non-equilibrium  plasmas  with  high  energy 
density.  The  software  is  being  designed  using  Object-Oriented  methods  for  maximal  flexibility,  portability 
and  easy  maintenance.  The  paper  describes  the  general  approach  used  in  the  model  development,  some  of 
the  problems  to  be  solved,  architecture  and  parallelization  strategies,  and  some  of  the  methods  currently  or 
being  implemented.  While  still  a  work  in  progress,  the  software  is  rapidly  building  capabilities  and  will  be 
applied  in  fundamental  simulations  of  various  plasma  discharges. 

1.  Introduction 

The  Air  Force  Research  Laboratory  has  a  significant  interest  in  the  development  and  evaluation  of 
transient,  high-power  and  micro-plasma  discharges  for  compact  high-performance  electric  propulsion 
systems,  high-power  conditioning,  novel  diagnostics  methods,  material  fabrication,  etc.  The  plasma 
conditions  in  these  new  generations  of  devices  are  usually  extremely  difficult  to  fully  characterize 
experimentally  or  theoretically.  Conventional  numerical  models  are  inadequate  in  those  regimes,  requiring 
the  use  of  “hybrid  models”  that  combine  different  numerical  models  for  different  components  of  the 

plasma,  along  with  the  interaction  with  electro-magnetic  fields  and  materials.  Hybrid  modeling  can  be 

particularly  challenging,  as  the  model  selection  can  be  dynamically  linked  to  the  plasma  conditions. 
Although  there  has  been  a  considerable  amount  of  research  in  developing  advanced  numerical  methods  and 

specialized  plasma  modeling  tools  for  specific  plasma  problems,  the  lack  of  an  integrated  suite  of  software 
tools,  combining  various  numerical  approaches,  is  preventing  the  investigation  of  some  types  of  problems. 
In  a  well-designed  software  environment,  one  should  be  able  to  use  the  same  software  tools  and  interfaces, 
irrespective  of  the  plasma  conditions  being  studied  (e.g.  grid  generation,  database,  transport,  visualization, 
diagnostics,  etc),  and  be  able  to  select  the  appropriate  methods,  or  allow  the  software  make  the  selection 
automatically.  The  current  research  effort  is  based  on  that  concept,  and  is  aimed  at  the  modeling  of  plasma 
discharges  characterized  by:  (a)  extreme  non-equilibrium  conditions;  (b)  small  spatial  and  temporal  scales, 
and;  (c)  high-energy  processes.  In  all  the  problems  of  interest  here,  the  methods  must  include  the 
interaction  with  local  and  applied  electromagnetic  (EM)  fields  and  collisional  momentum  and  energy 
exchange,  as  well  as  inelastic  and  radiative  processes. 

The  general  problem  reduces  to  the  evaluation  of  the  velocity  distribution  function  (DF)  that  describes  the 
system  and  its  dynamics,  i.e.  solving  the  Boltzmann  equation: 

2/  +  f  §j=cR(f,n  (i) 

ot 

where  f(t,x,  v)  is  the  velocity  DF  and  the  right-hand-side  (RHS)  symbolically  describes  all  exchange 
terms  due  to  collisional  or  radiative  effects.  However,  a  complete  knowledge  of  the  system  DF  is  not 
possible  given  the  current  and  even  foreseeable  computing  capabilities,  due  to  the  hyper-dimensionality  of 
the  phase  space,  especially  when  plasma  components  with  internal  structure  must  be  considered.  A 
practical  M&S  capability  for  the  multi-scale  problems  of  interest  here  can  only  be  achieved  through  hybrid 
modeling.  In  this  approach,  each  plasma  component  is  being  treated  by  the  numerical  method  that  is  the 
most  appropriate  and  most  efficient,  given  the  overall  physical  conditions. 
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The  plasma  discharges  of  particular  interest  here  are  all  characterized  by  a  high-energy  component  far  from 
equilibrium,  but  range  from  simple  planar  diodes  to  hollow  cathode  discharges,  vacuum  arcs,  and  intense 
laser-plasma  interactions.  The  planar  diode  is  an  excellent  verification  &  validation  case,  since  it  is 
extremely  simple  geometrically  (ID  problem),  yet  can  lead  to  a  rich  spectrum  of  dynamical  effects  which 
can  be  investigated  with  analytical  and  other  (Lagrangian)  methods  [1-2].  Another  excellent  test  bed  for 
hybrid  models  is  the  pseudospark  discharge,  which  has  been  previously  studied  as  a  high-power  plasma 
switch,  and  for  which  a  similar  attempt  at  analytical  description  of  the  ignition  phase  has  been  made  [3]. 

We  envision  the  integration  of  three  non-exclusive  categories  of  numerical  approaches:  (a)  continuum 
representations,  based  on  the  solution  of  partial  differential  equations  modeling  a  set  of  conservation  laws; 
(b)  particle  methods,  describing  the  plasma  as  a  collection  of  individual  “pseudo-particles”,  and;  (c)  direct 
kinetic  equation  solvers,  which  solve  the  dynamical  equations  (e.g.  Vlasov,  Fokker-Planck  and  Boltzmann 
equations)  in  a  high-dimensionality  phase  space.  The  latter  will  be  mostly  used  as  a  means  of  verifying  the 
hybrid  fluid-particle  model.  The  particle  models  include  both  Particle-In-Cell  (PIC)  and  Direct  Simulation 
Monte-Carlo  (DSMC).  The  PIC  model  [4]  is  initially  electrostatic,  non-relativistic  only,  and  will  be 
gradually  upgraded  to  a  fully  relativistic,  electro-magnetic  PIC.  The  coupling  between  fluid  and  stochastic 
solvers  proceeds  through  collisional  exchanges  of  mass  (e.g.  ionization),  momentum  and  energy;  these  are 
treated  through  a  Monte-Carlo  Collisions  (MCC)  algorithm,  and  variations  of  the  basic  procedure  are  being 
investigated.  The  DSMC  [5]  model  will  be  added  later.  The  fluid  models  currently  being  implemented  span 
a  range  of  non-equilibrium  descriptions,  starting  from  a  single-fluid  and  single-temperature  (thermal 
equilibrium)  model  all  the  way  to  multi-fluid  and  Collisional-Radiative  (C-R)  models. 

The  paper  is  organized  as  follows.  In  section  2,  the  approach  to  the  software  development  is  described, 
from  general  considerations  to  numerical  and  data  representation  issues,  language,  platform,  parallelization 
and  performance.  Section  3  describes  some  of  the  models  and  algorithms  implemented  to  date,  as  well  as 
some  integration  issues.  Section  4  describes  planned  future  work  and  summarizes  the  development  status. 

2.  Approach 

2.1.  General  Development  Strategy 

In  most  cases  of  interest  here,  the  plasma  can  be  decomposed  into  a  “bulk”  plasma  and  a  highly  energetic 
component  (electrons  and/or  ions).  If  the  bulk  plasma  is  strongly  collisional  (i.e.  sufficiently  dense),  it  can 
be  modeled  as  a  fluid,  albeit  at  various  levels  of  thermal  non-equilibrium.  The  energetic  component  can 
itself  be  a  function  of  time  and  space;  for  example,  at  the  onset  of  the  discharge  there  can  be  very  few 
energetic  particles;  during  a  mid-phase  the  energetic  component  is  maximum,  until  it  finally  relaxes 
through  collisional  exchange  with  the  bulk  and  disappears.  One  could  also  envision  a  reverse  process, 
where  energetic  particles  are  generated  from  the  bulk  (e.g.  stochastic  heating)  itself.  Short  of  modeling  the 
entire  DF  in  phase  space,  one  would  need  to  devise  methods  for  transferring  information  back  and  forth 
between  the  bulk  and  energetic  components.  This  is  the  general  problem  of  mass,  momentum  and  energy 
coupling  between  the  dual  representations  (fluid  and  particles). 

This  issue  is  exemplified  in  previous  attempts  at  pseudospark  modeling  with  a  hybrid  method  [6];  in  that 
model,  an  energy  threshold  is  assumed  (the  ionization  potential)  such  that  the  plasma  is  modeled  as  a  fluid 
below  it,  and  as  particles  above  it.  However,  there  are  some  potential  problems  with  this  approach.  First, 
there  is  no  definite  energy  threshold  of  behavior:  if  the  threshold  is  the  ionization  potential,  this  implies  that 
all  electronic  excitations  are  being  ignored.  If  it  includes  the  minimal  gap  for  electronic  excitation,  this  does 
not  include  excitation  of  other  modes  for  molecular  plasmas.  Furthermore,  that  particular  model  did  not 
treat  new  generations  of  electrons  (liberated  by  ionizing  collisions)  as  particles,  and  while  the  ballistic 
approximation  is  efficient  in  separating  the  time  scales  of  the  electrons  and  ions,  it  does  not  allow  for 
immediate  field-coupling  effects  on  the  high  energy  electrons,  which  might  be  necessary  for  the  beam 
phase.  It  is  clear  that  a  better  hybrid  model  may  be  necessary  in  order  to  successfully  model  all  phases  of 
the  pseudospark  operation,  as  well  as  other  high-energy  and/or  micro-discharges,  and  this  may  require  that 
a  consistent  mechanism  for  particle  creation  and  absorption  from  and  into  the  bulk  be  devised. 

Another  problem  to  address  concerns  the  collisional  models.  One  can  model  collisions  between  the 
energetic  particles  and  the  background  gas  through  a  Monte  Carlo  Collisions  (MCC)  procedure.  For 
electron  collisions,  the  procedure  is  similar  to  a  test-particle  Monte-Carlo  [7]:  due  to  the  difference  in 
masses,  the  target  atom  can  be  assumed  unperturbed  during  the  collision.  In  the  case  where  the  energetic 
particle  is  an  ion  itself,  particles  must  be  sampled  from  the  equilibrium  bulk  distribution,  collided  with  the 
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energetic  particles,  and  the  momentum  and  energy  exchanged  monitored;  this  approach  does  exactly 
conserve  momentum  and  energy,  but  is  expensive.  Another  approach  is  to  use  the  more  efficient  “collision- 
field”  method  [8-10],  which  can  also  conserve  momentum  and  energy  but  averages  the  velocity  dependence 
into  a  collision  frequency;  this  approach  can  also  work  for  electron  collisions.  Finally,  there  is  the 
interesting  alternative  of  extended  pseudo-particles,  each  one  carrying  its  own  internal  distribution  at  an 
individual  “temperature”  [11].  The  standard  MCC  is  currently  being  implemented;  at  a  later  date,  we  plan 
to  also  implement  the  collision-field  approach  and  compare  the  results. 

It  is  also  important  to  be  able  to  validate  the  hybridization  and  collisional  methods  implemented  in  the 
proposed  code;  for  that,  we  will  also  develop  and  implement  some  simple  Vlasov  and  Fokker-Planck 
kinetic  solvers  with  high  resolution  in  phase  space.  This  will  allow  us  to  test  procedures  for  the  exchange  of 
particle  and  bulk  components  (without  collisions).  Furthermore,  we  will  also  implement  a  multi-fluid 
solver,  which  will  model  the  electrons,  ions  and  neutrals  separately  as  fluids,  with  friction  and 
electromagnetic  coupling.  This  will  allow  us  to  test  the  limits  at  which  the  fluid  model  starts  to  break  down. 

Given  the  disparity  of  the  modeling  approaches  and  the  complexity  of  the  physics,  it  does  not  appear 
possible  that  all  approaches  can  be  combined  into  a  single  code.  Such  a  code  would  be  extremely  large; 
very  difficult  to  use  even  on  simple  problems;  extremely  difficult  to  maintain  and  upgrade,  and  using 
massive  amounts  of  computational  resources.  It  is  much  more  appropriate  to  develop  a  library  of  modeling 
tools,  specialized  for  plasma  simulations,  along  with  the  software  required  for  the  modeling  of  the  specific 
applications  of  interest.  Therefore,  other  researchers  can  import  these  tools  or  integrate  them  to  construct 
another  code  designed  for  another  set  of  applications.  In  addition,  the  software  should  have  a  number  of 
features  aimed  at  maximizing  its  usefulness,  such  as:  multi-dimensional  capability  (ID,  2D,  3D);  ability  to 
handle  complex  shapes  and  boundary  conditions;  simplified  maintenance  and  upgrade;  portability;  high 
computational  efficiency,  and;  graphical  user  interface  (GUI)  and  visualization  tools.  These  objectives  are  a 
part  of  the  overall  software  development  strategy. 

2.2.  Software  Architecture  and  Parallelization  Strategy 

To  have  a  significant  impact  on  problems  of  practical  interest,  one  must  be  able  to  efficiently  combine  the 
most  important  types  of  numerical  approaches  used  for  plasma  modeling.  Therefore,  one  must  construct  the 
data-structures  and  associated  methods  to  be  used  in  the  software  framework:  these  include  geometry,  flow, 
field  variables,  and  particles.  Generalized  classes  of  solvers,  easily  re-configurable,  must  be  constructed, 
and  the  coupling  physics  and  interface  procedures  must  be  implemented.  For  a  number  of  reasons,  Java  has 
been  chosen  as  the  language  for  the  design  and  development  of  the  software  architecture.  Java  is  a  strong 
Object  Oriented  (OO)  design  language  with  an  easy  syntax  that  allows  us  to  easily  develop,  debug  and 
share  software  modules.  Portable  visualization  and  GUI  can  be  completely  integrated  using  open-source 
software,  and  there  are  no  issues  with  compiler  compatibilities  between  different  vendors  and  platforms. 
Furthermore,  Java  allows  us  to  easily  test  various  parallelization  approaches  across  different  platforms,  in  a 
way  that  is  compatible  with  the  modem  approach  of  “grid-computing”.  Using  a  combination  of  Remote 
Method  Invocation  (RMI)  and  multi-threading,  one  can  examine  both  distributed  and  shared  memory 
parallelization,  and  have  remote  steering  and  real-time  visualization.  Except  for  a  few  exceptions  [12-14], 
Java  has  so  far  not  been  considered  as  a  viable  alternative  for  large-scale  scientific  computing  because  of 
the  perceived  lack  of  performance.  However,  recent  versions  of  the  Java  Virtual  Machine  (JVM)  have 
come  a  long  way  to  improve  this  situation;  our  own  tests,  confirmed  by  third-party  developers,  indicate  that 
at  least  on  some  platforms  (Intel)  the  execution  speed  for  simple  floating-point  calculations  are  similar  to 
optimized  C.  More  complete  performance  comparisons  have  been  performed  elsewhere  with  similar 
conclusions  [15].  Differences  of  even  10-30%  in  execution  speed  may  not  be  a  critical  factor  when 
compared  to  other  issues  such  as  platform  availability  and  cost,  whether  hardware  (the  ability  to  use 
existing  networks  of  disparate  workstations  has  a  significant  cost  advantage,  especially  given  the  rapid 
obsolescence  of  computer  technology)  and  personnel  costs  for  development  and  maintenance.  A  later 
version  of  the  software  may  be  written  in  C++  if  it  turns  out  to  be  beneficial  -  the  translation  from  Java  to 
C++  is  not  expected  to  present  major  difficulties.  It  is  also  expected  that  the  simulations  will  be 
computationally  demanding  (especially  in  3D);  therefore,  it  is  necessary  to  implement  parallelization 
strategies  as  early  as  possible  (once  the  data-structures  and  code  architecture  have  been  stabilized),  while 
maintaining  a  high  degree  of  flexibility. 
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There  are  several  levels  of  parallelization  that  could  be  applied  to  the  code.  The  current  trends  in  computing 
and  communication  technologies  increasingly  favor  the  development  of  clusters  of  computing  nodes,  each 
node  having  an  increasing  number  of  shared-memory,  hyper-threaded  processors.  Such  clusters  are 
routinely  available  with  up  to  128  nodes  of  2  CPU  each,  with  specialized  machines  having  thousands  of 
nodes.  Multi-processor  workstations  are  also  readily  available  platforms  that  can  be  tapped  for  low-cost 
computing.  It  is  therefore  likely  that  an  efficient  parallelization  strategy  should  plan  for  a  re-configurable 
network  of  multiple  platforms  and  a  hierarchy  of  distributed  and  shared-memory  architectures;  this  has 
been  a  major  factor  in  choosing  Java  as  the  top-level  language  for  the  software  architecture.  Drawing  from 
the  “grid-computing”  paradigm,  we  have  created  a  client-master-server  architecture  (see  Figure  2)  such  that 
computations  can  be  started  on  any  available  platforms,  and  monitored  in  real-time  by  the  user/client.  This 
is  similar  to  the  JavaGrid  approach  described  in  [12],  albeit  at  a  lower  level  of  complexity  and  capability. 

To  keep  maximum  flexibility  in  mind,  the  code  is  designed  to  make  inter-node  communication  as 
transparent  as  possible.  This  is  accomplished  by  the  implementation  of  Java’s  remote  method  invocation 
(RMI)  for  all  inter-processor  communications  [16].  Java-RMI  allows  for  the  manipulation  of  remote 
objects  on  a  local  machine.  As  an  example,  a  client  can  invoke  remote  methods  on  the  master,  asking  the 
remote  method  to  return  an  object  to  the  client,  such  as  the  computation  data  for  visualization.  The  RMI 
architecture  of  the  code  consists  of  three  levels — client,  master,  and  server.  The  code  is  initialized  on  the 
master  with  a  list  of  server(s)  that  will  handle  the  computations.  A  “server”  is  a  single-term  denomination 
that  can  physically  be,  for  example,  a  single-processor  or  multi-processor  workstation,  or  a  node  in  a  Linux 
cluster*.  Once  the  master  has  been  initialized,  the  client  can  establish  a  connection  to  it  \  The  client  can  then 
perform  preprocessing  operations  such  as  set-up  of  the  grid,  boundary  conditions,  chemistry,  etc,  to 
initialize  the  computations  and  engage  in  real-time  visualization  of  the  data  as  it  is  being  computed. 

The  software  of  course  supports  configurations  from  single-processor  workstations  to  multi-processor 
shared-memory  workstations  as  well  as  multi-processor  distributed-memory  clusters,  in  compliance  with 
the  grid-computing  paradigm.  This  is  made  possible  by  the  multi-threading  capability  of  the  server  code, 
which  allows  for  multiple  blocks  (each  block  being  a  subset  of  the  computational  domain)  to  be  computed 
on  each  “server”.  Each  block  runs  in  its  own  thread,  and  as  a  result,  multiple  blocks  on  a  single  processor 
are  computed  in  parallel.  Multiple  servers  can  be  initialized  with  multiple  blocks  residing  on  each  server. 
Only  one  server  instance  is  needed  per  computer  node,  and  as  many  servers  can  be  initialized  as  there  are 
blocks.  The  computation  iterations  are  carried  out  on  each  server  along  with  inter-server  communications. 

The  synchronization  of  all  iterations  and  communications  is  the  responsibility  of  the  master.  The  master 
utilizes  various  locks  and  the  associated  wait/notify  methods  inherited  from  java .  lang.Ob  ject  to 
determine  when  all  of  the  servers  have  finished  their  computation  and  communication  iterations  and  to  cue 
them  so  that  they  may  proceed  with  the  next  iteration.  Once  all  of  the  blocks  have  been  locked,  the  master 
is  notified,  and  in  turn,  notifies  all  the  blocks  to  begin  the  next  iteration.  The  master  is  also  in  charge  of 
global  data  which  is  processed  from  local  data  from  each  block,  such  as  the  time  step,  or  data  that  must  be 
visualized.  The  client  code  provides  a  graphical  user  interface  from  which  the  code  can  be  controlled  and 
data  can  be  viewed  in  real-time.  The  client  can  make  requests  to  the  master  to  grab  data  from  the  servers  as 
it  is  being  computed  and  display  it  in  real-time  with  VisAD — a  Java  visualization  API  built  upon  Java3D 
[17].  It  can  be  selected  from  the  client  as  to  which  blocks  to  view  data  from  so  that  the  communication 
overhead  does  not  bog  down  the  overall  computation  time. 

The  most  crucial  concern  regarding  communication  overhead  is  the  optimization  of  the  intra-/inter-server 
communications.  At  every  time- step,  it  is  necessary  for  each  block  in  the  domain  to  share  its  information 
with  one  or  more  neighboring  blocks.  If  two  neighboring  blocks  are  being  computed  on  the  same  server, 
the  communication  can  be  handled  by  standard  local  method  calls.  However,  if  two  neighboring  blocks  are 
computed  on  different  servers,  then  inter-server  communication  is  involved  and  it  becomes  necessary  to  use 
Java-RMI  or  other  suitable  message  passing  protocol.  Figure  3  illustrates  how  two  neighboring  blocks 
communicate  with  each  other.  There  are  two  types  of  patches  which  allow  data  transfer  between  two 


’’  One  should  think  of  the  “server”  as  a  computational  unit  that  provides  a  “service”  (computations).  There  is  a-priori  no 
reason  why  the  “server”  cannot  by  itself  further  distribute  the  computational  load  to  other  service  providers,  and 
although  this  recurrent  (multi-level)  parallelization  is  not  yet  implemented,  this  is  the  direction  taken  here. 

*  The  master  is  a  command  and  control  center,  responsible  for  gathering  and  distributing  information  to  the  servers  and 
client,  scheduling  the  tasks,  providing  critical  I/O,  and  monitoring  commands  from  the  client. 
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blocks.  The  first  of  these  being  the  source  patch  is  an  array  of  objects  which  contain  references  to  data 
available  in  the  boundary  layer  of  cells.  The  data  references  are  depicted  by  the  dashed  arrows  in  Figure  3. 
The  second  is  the  target  patch  from  which  the  ghost  layer  of  cells  receives  its  information.  During  block 
communication,  there  is  a  method  call  for  each  source/target  patch  pair  in  which  the  target  patch  is  filled 
with  the  information  from  its  neighboring  block’s  source  patch.  The  method  calls  are  depicted  by  the  solid 
arrows  which  can  either  be  a  local  call  if  both  blocks  are  on  the  same  server  or  an  RMI  call  if  each  block  is 
on  a  separate  server.  The  exchange  of  boundary  condition  information  between  servers  is  done  without 
involving  the  master:  this  prevents  a  bottle-neck  when  scaling  to  a  large  number  of  processors. 

2.3.  Geometry  Data-Structure 

The  first  compatibility  issue  one  must  deal  with  concerns  the  grid  structure.  Ideally,  the  software  should  be 
grid-independent;  this  implies  that  one  could  at  first  decide  to  implement  solvers  based  on  structured  grids, 
while  reserving  the  ability  to  implement  unstructured  grids  and  their  solvers  at  a  later  stage.  All  numerical 
approaches  described  above  make  use  of  a  discretized  space:  besides  the  self-evident  case  of 
Computational  Fluid  Dynamics  (CFD),  PIC  is  a  typical  particle-mesh  method,  and  DSMC  needs  a  mesh  to 
determine  collision  partners.  Grids  can  be  described  with  the  following  categories:  (1)  Structured  (constant 
number  of  grid  points  in  each  direction),  or  (2)  Unstructured  (triangles,  tetrahedra).  The  former  also  has  the 
following  sub-cases:  (1-a)  Cartesian  (Ax,  Ay,  Az  are  constants);  (1-b)  Uniform  (spacings  are  not  constants); 
(1-c)  Curvi-linear,  or  “body-fitted”  (not  necessarily  orthogonal).  Structured  grids  can  be  described  by 
multi-dimensional  arrays  of  given  dimensions  (Q  [  ]  [  ]  [  ]  ),  i.e.  any  grid  point  can  be  addressed  by  a  set  of 
indices  (i,j,k  in  3D).  In  unstructured  grids  there  is  no  such  ordering,  and  the  grid  connectivity  (i.e.  how  to 
access  another  point  from  the  current  one)  requires  a  list  of  pointers  between  grid  elements.  Despite  the 
apparent  difference  between  the  two  descriptions,  one  can  use  the  same  approach  in  representing  the  grid 
data  and  connectivity.  Thus,  one  can  define  an  object  gNode  to  represent  a  grid  point  (a  node  in  the  mesh); 
each  gNode  contains  the  coordinates  of  the  point  (in  real  or  transformed  space),  as  well  as  a  list  of  pointers 
(addresses  in  the  memory  where  the  data  is  located,  or  integer  indices  that  indicate  a  position  in  an  array)  to 
the  neighboring  points.  The  same  goes  for  object  of  type  “  gCeli’,  which  describes  the  volume  enclosed  by 
the  grid  lines  or  grid  surfaces.  Thus  one  can  walk  through  the  connected  list  of  nodes  and  cells  to  reach 
every  point,  whether  the  grid  is  structured  or  not.  Of  course,  the  structured  case  has  the  additional 
advantage  that  one  could  also  directly  address  a  given  node  or  cell  from  the  set  of  indices  ( i,j,k).  However, 
instead  of  forcing  the  data  to  be  stored  in  a  structured  form,  it  is  just  as  easy  to  have  a  method  that  provides 
the  memory  location  of  a  node/cell  with  those  particular  indices.  Thus,  the  data  structure  (a  one¬ 
dimensional  list  of  objects)  can  remain  the  same  to  describe  both  structured  and  unstructured  grids,  while 
the  additional  methods  “ getNode ”  and  “ getCell ”  allow  a  mapping  between  the  structured  indices  (i,j,k), 
arguments  to  the  method,  and  the  list.  Such  a  mapping  is  of  course  void  in  the  unstructured  case. 

Within  this  approach  a  “Mesh”  object  is  defined  as  a  list  (or  ID  array)  of  cells,  nodes  and  surfaces,  with 
pointers  to  connected  objects.  If  a  structured  solver  is  being  used,  one  must  address  information  on  either 
side  of  a  cell  or  surface  by  the  pointers  to  the  corresponding  objects,  e.g.  in  2D,  U  [  i-1  ]  [  j  ]  is  obtained 
from  U[ncL]  ,  with  ncL=pointer_left_cell  (nc).  This  procedure  has  two  disadvantages  with 
respect  to  a  simple  multi-dimensional  data  structure  (e.g.  U  [  ]  [  ] ):  (a)  there  is  redundant  information  being 
stored  in  the  Mesh  object;  however,  the  information  is  contained  in  a  few  variables  and  does  not  constitute 
any  significant  penalty;  (b)  there  is  additional  overhead  in,  for  example,  the  method  call  “getCell”,  which 
maps  the  one-dimensional  index  of  the  list  of  cells  with  the  structured  indices  (i,j,k)',  however,  one  can 
easily  construct  algorithms  such  that  this  method  call  is  relatively  infrequent  and  the  cost  is  minimal' 
compared  to  the  large  number  of  operations  involved  in  a  flow  or  field  solver,  for  example.  The  major 
advantage  is  that  this  approach  allows  the  same  data- structure  to  be  used  for  any  type  of  grid  structure.  The 
Mesh  object  can  be  sufficiently  general  that  each  cell  has  an  arbitrary  connectivity,  while  some  variables 
used  only  for  some  types  of  grids  may  remain  un-instantiated  for  others.  For  example,  a  structured  grid  may 
need  a  metric  for  some  solvers;  this  metric  may  be  computed  and  stored  in  that  case,  but  is  empty  when  the 
grid  is  unstructured.  Similarly,  some  techniques  of  dynamic  grid  adaptation  may  be  valid  only  for  some 
types  of  grids  (e.g.  octree  subdivision  for  Cartesian  grids)  and  not  on  others.  Each  Mesh  object  is  therefore 
flagged  by  a  few  Boolean  variables  (e.g.  Cartesian?  Curvilinear?  Unstructured?),  so  that  one  can  match 
appropriate  solvers  and  methods  to  the  types  of  grids  being  used.  Even  if  not  all  solvers  are  immediately 


1  In  some  languages,  in-lining  can  be  used  to  eliminate  the  overhead  associated  with  a  function  call. 
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implemented,  future  extensions  are  facilitated  because  the  geometry  infrastructure  can  support  it.  This  is 
the  approach  used  here:  while  the  current  focus  is  on  developing  solvers  for  structured,  body-fitted  grids, 
unstructured  grid  implementations  will  be  facilitated  in  the  future,  using  solvers  and  grid  generators  (which 
must  provide  the  connectivity  information)  from  external  sources. 

Structured,  body-fitted  grids  are  efficient  for  CFD,  especially  when  boundary  layers  are  present;  in  that 
case,  cells  can  be  stretched  to  match  the  preferred  direction  of  the  gradients  (normal  to  the  body  surface), 
while  allowing  efficient  tri-diagonal  implicit  solvers  and  increasing  the  accuracy  of  dimensionally-split 
high-order  schemes.  This  approach  may  also  be  efficient  for  plasma  sheaths,  since  again  the  thermal, 
plasma  composition  and  plasma  potential  gradients  are  aligned  with  the  normal  to  the  body  surface.  On  the 
other  hand,  Cartesian  grids  are  preferable  for  particle  transport,  the  construction  of  simple  high-order 
schemes  (e.g.  for  the  Maxwell  equations),  Fourier-transformed  schemes,  and  dynamic  cell  adaptation.  If  the 
grid  is  Cartesian,  it  is  extremely  easy  to  keep  track  of  the  cell  to  which  the  particle  belongs,  given  its 
position,  e.g:  i=x/Ax.  This  can  be  important  for  a  large  number  of  non-interacting  particles,  where 
transport  (particle  “pusher”)  is  the  major  computational  activity.  However,  this  is  true  mostly  for  collision¬ 
less  plasmas  in  simple  geometries1)  When  the  plasma  is  collisional,  most  of  the  computational  cost  is  in  the 
MCC  procedure.  In  those  cases,  the  additional  cost  in  determining  the  relative  particle  position  with  respect 
to  the  cells  (i.e.  cell-indexing  the  particle)  for  non-Cartesian  geometries  becomes  a  negligible  fraction.  This 
extra-cost  is  typically  a  few  arithmetic  operations,  and  is  roughly  of  the  same  order  for  body-fitted  grids 
(using  metric-transformed  variables)  and  unstructured  grids  (computing  relative  distances),  as  long  as  the 
particle  moves  at  most  to  neighboring  cells.  This  last  limitation  is  also  not  a  restriction  for  collisional 
plasmas,  since  it  is  important  for  solution  accuracy.  One  should  point-out  that  the  two  approaches  can  be 
effectively  combined,  by  using  body-fitted  grids  near  the  body  and  in  collisional  regions,  and  Cartesian 
grids  away  from  the  body  and  in  collisionless  regions.  The  ability  to  handle  any  grid  type,  as  described 
above,  is  essential  to  this  optimization  scenario. 

2.4.  Fluid  Data-Structure 

The  data-structure  for  the  fluid  and  plasma  variables  must  also  be  highly  flexible.  The  fluid  can  be 
described  at  various  levels  of  non-equilibrium,  corresponding  to  various  relaxation  time-scales.  By 
assumption,  the  velocity  distribution  of  the  fluid-like  plasma  components  remains  Maxwellian,  the 
collisional  momentum  exchange  between  particles  of  similar  masses  being  the  fastest  relaxation  process. 
However,  electrons  and  ions  may  have  different  velocity  distributions  (i.e.  different  translational 
temperatures),  and  internal  modes  may  not  follow  identical  equilibrium  distributions,  or  may  significantly 
deviate  from  Boltzmann  distributions.  This  hierarchy  of  non-equilibrium  descriptions  results  in  a  variety  of 
multi-temperature  models.  The  hierarchy  for  atomic  plasmas  would  be  single-temperature,  two-temperature 
(Th  -Te ),  three-temperature  (T/*  -Te  -Tx),  and  Collisional-Radiative  (C-R)  models,  where  7/*  is  the 
translational  temperature  of  heavy  particles,  Te  is  the  temperature  of  the  free  electrons,  and  Tx  is  the 
temperature  of  the  (assumed)  Boltzmann  distribution  of  electronic  excited  states.  In  the  two-temperature 
model,  the  bound  electronic  states  are  assumed  to  be  in  a  Boltzmann  distribution  at  the  temperature  Te  ;  this 
is  usually  valid  in  cases  where  the  electron  collisions  dominate  over  atomic  collisions.  In  low-density  cases, 
when  radiative  rates  are  comparable  to  collisional  or  atomic  quenching  becomes  significant,  the  electronic 
states  cannot  be  described  by  a  Boltzmann  distribution;  in  this  regime  the  C-R  model  is  required. 

The  ability  to  generate  a  variety  of  non-equilibrium  descriptions  is  another  important  aspect  of  the  code 
flexibility.  One  can  determine  from  the  input  conditions  which  model  is  necessary,  and  also  easily  perform 
parametric  studies  to  verify  assumptions  regarding  the  state  of  the  plasma,  or  evaluate  the  importance  of 
non-equilibrium  effects.  Furthermore,  it  is  also  desirable  in  some  circumstances  to  adjust  the  level  of 
description  according  to  the  local  conditions.  This  would  maximize  efficiency,  since  higher  levels  of  non¬ 
equilibrium  imply  higher  computational  costs.  The  software  is  currently  being  configured  to  allow  different 
non-equilibrium  fluid  models  in  different  regions;  this  feature  is  currently  limited  to  static  and  pre¬ 
determined  conditions,  but  it  is  theoretically  feasible  to  extend  this  capability  to  automatic  and  time- 
dependent  determination  of  the  non-equilibrium  model.  This  capability  implies  that  boundary  conditions 
between  contiguous  domains  automatically  transform  the  data  between  lower  and  higher-levels  of 


When  complex  3D  shapes  are  involved,  the  Cartesian  cells  must  be  “cut”  near  the  boundary,  leading  to  complex 
tetrahedral  shapes.  To  achieve  high  accuracy  near  the  boundaries  becomes  more  difficult  and  efficiency  drops 
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thermodynamic  description.  Therefore,  key  methods  to  be  implemented  in  the  library  for  fluid  modeling  are 
for  the  contraction  (from  higher-  to  lower  level)  and  expansion  of  information.  An  example  of  the  former 
would  be  the  grouping  of  all  excited  levels  from  a  C-R  model  into  a  Boltzmann  distribution  and  the  forcing 
of  equilibrium  with  the  free  electrons  (from  C-R  to  2-T  model);  an  example  of  the  latter  would  be  the 
generation  of  an  additional  variable  for  the  electron  energy  from  the  single  fluid  equations  (from  1-T  to  2- 
T).  Preliminary  tests  of  these  procedures  have  been  performed  successfully.  Of  course,  the  procedures  are 
accurate  as  long  as  the  conditions  permit  it,  i.e.  are  such  that  the  higher  level  of  description  (e.g.  2-T  model) 
is  compatible  with  local  thermodynamic  equilibrium  (LTE)  for  the  mode  being  expanded  or  contracted 
(Te  ~  T/t ).  This  approach  yields  maximal  benefits  when  the  high  level  of  thermodynamic  description  is 
very  expensive  to  maintain  in  near-equilibrium  conditions;  this  would  be  the  case  of  the  C-R  model,  for 
example,  since  it  requires  a  solver  for  the  master  equations  of  collisional  and  radiative  excitation  and  de¬ 
excitation,  and  mass  conservation  equations  for  all  levels.  Similar  benefits  could  be  found  for  molecular 
plasmas.  Note  also  that  a  similar  procedure  applied  to  translational  modes  (i.e.  velocity  DF)  would  allow 
the  transition  between  CFD  in  dense  fluid  regions  and  DSMC  in  transitional  or  rarefied  regions. 

Finally,  one  should  also  point  out  that  atomic  electronic  states  are  too  numerous  to  be  treated  independently 
inside  a  C-R  model;  the  usual  procedure  consists  of  grouping  together  electronic  states  which  are  close 
together  in  the  energy  spectrum. 

2.5.  Particle  Data-Structure 

Particles  usually  have  a  very  simple  structure:  they  must  contain  a  minimum  of  information,  since  accurate 
simulations  requiring  106-107  particles  or  more  would  consume  an  inordinate  amount  of  memory  if  each 
particle  object  contained  an  excessive  amount  of  data.  The  minimum  information  is:  (a)  a  particle  type 
identifier  (this  can  be  the  species  number);  (b)  particle  coordinates  (x,y,z);  (c)  particle  velocity  (vx,vy,vz). 
The  rest,  i.e.  charge,  mass,  etc,  can  be  obtained  from  the  species  information.  Additional  information  may 
include:  (d)  particle  statistical  weight;  (e)  cell  number  to  which  the  particle  belongs,  but  also;  (f) 
transformed  local  coordinates  (£,r|,£)  which  are  especially  useful  for  monitoring  whether  a  particle 
crosses  into  an  adjacent  cell  (and  which  one)  in  body-fitted  grids,  and;  (g)  time  scale  -  this  will  be 
explained  further  in  section  3.1.  The  statistical  weight  (each  simulated  particle  is  a  “pseudo-particle”  and 
represents  a  large  number  05  of  real  particles)  could  be  made  species-dependent  only,  in  which  case  it 
would  not  need  to  be  included  as  item  (d)  in  the  particle  data- structure,  but  there  are  several  reasons  why 
this  may  be  beneficial:  (1)  when  particles  are  injected  from  time-dependent  boundary  conditions  (e.g. 
electrons  from  cathodes)  a  higher  accuracy  could  be  obtained  by  varying  the  statistical  weight  in  addition 
to,  or  instead  of,  varying  the  number  of  pseudo-particles;  (2)  chemical  reactions  can  be  more  accurately 
computed  using  variable  statistical  weights;  (3)  higher  statistical  accuracy  can  be  achieved  when  the 
particle  density  varies  by  several  orders  of  magnitude  across  the  computational  domain,  by  duplicating 
particles  that  move  into  low-density  regions.  These  issues  are  common  to  DSMC  simulations,  where  the 
use  of  variable  statistical  weights  is  very  frequent  [5]. 

Since  some  of  these  quasi-particles  may  describe  gas  species  with  internal  structure  (electronic,  vibrational 
and  rotational  states),  there  may  also  be  some  information  required  to  describe  these  modes.  Ideally,  to  each 
particle  would  be  associated  a  unique  state;  however,  to  accurately  (i.e.  with  low  statistical  noise)  describe 
the  complete  DF,  one  would  require  an  extraordinarily  large  number  of  particles.  One  must  instead  attribute 
a  distribution  of  internal  states,  e.g.  a  simulated  molecule  could  have  its  own  rotational  temperature, 
indicating  how  these  modes  are  populated.  During  a  collisional  energy  exchange,  one  can  compute  and 
modify  the  overall  internal  energy  associated  with  that  mode,  and  assign  again  a  temperature.  The 
disadvantage  of  this  approach  is  that  one  potentially  neglects  correlations  between  internal  and  translational 
modes.  For  example,  if  situations  occur  where  most  of  the  excited  states  (indicated  by  a  quantum  number 
A )  are  associated  with  particles  moving  at  high-velocity  in  the  z-direction,  the  statistical  correlation 

(A|  vz  )  can  be  high  when  both  A  and  Vz  are  high,  and  this  effect  would  not  be  reproduced  if  internal 
energy  distributions  are  stored  for  each  particle.  Such  situations  can  occur  in  special  cases  (beam 
interactions)  and  although  molecular  plasmas  are  not  currently  being  considered,  it  is  prudent  to  design  the 
software  with  both  potential  approaches  in  mind;  the  data- structure  can  simply  have  un-instantiated  arrays 
of  levels  or  distribution  parameters,  which  can  be  used  for  whichever  model  becomes  implemented. 
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This  approach  also  could  be  extended  to  the  velocity  distribution,  i.e.  the  particle  data  includes  also  a 
velocity  distribution.  This  approach  is  particularly  interesting  when  dealing  with  the  problem  of  particle 
coalescence,  to  be  described  in  the  next  section. 


3.  Numerical  Model  Development 
3.1  PIC  Model 

The  PIC  model  currently  in  place  is  a  simple  electrostatic,  non-relativistic  model.  Particles  are  generated 
with  variable  statistical  weight;  this  is  useful  in  cases  where  emission  is  variable,  i.e.  depends  on  the  plasma 
and  field  conditions  near  the  cathode.  Only  electrons  are  presently  implemented  in  the  PIC  model,  and 
electrode  boundary  conditions  are  limited  to  a  couple  of  simplified  models.  Electrons  are  “pushed”  in  three 
partial  steps:  first,  the  velocity  is  advanced  to  the  time  level  (n+1/2);  second,  the  particle  position  is 
changed  using  the  updated  velocity;  third,  the  particle  velocity  is  advanced  again  to  the  next  time  level 
(n+1).  This  is  the  simplest  form  of  a  second-order  symplectic  integration  scheme  [18].  Note  that  the 
velocity  must  be  evaluated  at  time  level  (n+1)  if  one  wants  to  evaluate  the  current  density  and  inelastic 
collision  rates  appropriately  at  that  level.  The  global  time  step  is  of  the  order  of  the  average  time  taken  by  a 
particle  to  move  by  the  characteristic  length  of  a  cell.  There  is  no  strict  Courant-like  condition  for  the 
particles,  but  that  in  average  the  particles  will  not  cross  more  than  one  cell  per  time-step.  Exception  is  made 
for  the  particles  in  cells  lying  next  to  a  boundary,  as  explained  below. 

The  “move”  phase  (second  step  of  the  list  above)  is  in  itself  decomposed  into  several  phases.  As  noted  in 
section  2.5,  each  particle  data-structure  contains  a  time  variable  (item  (g)),  which  is  the  time  left  to  the 
particle  for  executing  its  move.  At  the  start  of  each  iteration,  this  time-left  is  increased  by  the  global  time 
step  At .  During  the  first  phase,  only  the  particles  that  belong  to  cells  near  a  boundary  are  being  moved 
(this  is  a  cell-based  operation,  which  obviously  requires  that  the  cell  contains  a  list  of  particles  inside  the 
cell).  For  these,  the  exact  time  at  which  the  particle  hits  the  boundary  (if  at  all)  is  determined,  and  these 
particle  positions  are  advanced  up  to  that  time;  this  implies  that  the  time-left  is  decreased  by  the 
corresponding  amount.  These  particles  are  then  flagged  and  treated  by  the  particle  boundary  condition 
procedure.  For  example,  if  the  boundary  is  a  solid  wall,  the  particle  velocity  is  reflected,  and  if  the 
boundary  is  a  patch  to  another  grid,  the  particle  is  buffered  for  transmission  to  another  processor.  Finally, 
the  particles  are  moved  again  by  the  remaining  time-left  individual  to  each  particle.  This  last  phase  includes 
a  check  on  the  new  cell  position;  if  it  appears  that  some  particle  crosses  into  a  boundary  (i.e.  was  two-layers 
removed  from  the  boundary,  and  therefore  not  considered  by  the  cell-based  pre-positioning),  its  trajectory 
is  re-computed  with  a  smaller  elapsed  time  until  it  end-up  within  the  computational  domain.  This  procedure 
prevents  the  loss  of  particles  from  the  system;  it  would  affect  only  a  few  particles,  and  only  if  the  global 
time  step  is  chosen  relatively  large,  and  therefore  is  not  a  significant  computational  penalty.  Note  that  such 
particles  would  end-up  with  some  time  remaining,  which  will  be  processed  at  the  next  time-step.  The  use  of 
a  time-left  variable  for  each  particle  allows  some  flexibility  in  the  particle  transport,  avoiding  the 
processors  to  spend  extraneous  time  attempting  to  synchronize  all  particles.  Note  also  that  the  check  on  cell 
identification  from  the  particle  position,  at  the  end  of  the  final  “move”  phase,  is  a  Mesh- based  method  that 
is  dependent  on  the  grid  structure.  In  the  Cartesian  case,  a  few  simple  divisions  provide  the  answer;  in  the 
body-fitted  case  or  unstructured  grid  case,  additional  operations  may  be  involved. 


A  critical  issue  in  the  problems  of  interest  here  is  the  high  rate  of  ionization  that  can  be  obtained.  If  each 
ionization  event  from  a  simulated  collision  of  a  particle  with  atoms  of  the  background  fluid  generates 
additional  particles,  the  number  of  simulated  particles  can  grow  exponentially,  leading  to  a  commensurate 
growth  in  computational  work.  To  prevent  this,  one  must  devise  an  algorithm  for  particle  coalescence.  Fet 
us  first  consider  the  merging  of  two  particles  of  identical  type  (and  mass  m)  into  a  single  pseudo-particle  of 
combined  statistical  weight  03a£=G3a4€iJ£  ,  and  with  a  combined  velocity  uab=$5  a  Pa  ^  a  Pa  )  /(mU5  ab ) . 
This  procedure  conserves  momentum,  but  the  sum  of  the  kinetic  energies  of  the  individual  constituents 
differs  from  the  mean  kinetic  energy: 

p  k 

SE~EKa  +EKb~EKab  ~m a  +L+t5fe  ~  ~m  ab  ~ 

2m  2m  2 


ft2 

Uab 


Note  that  8F  is  always  positive: 


®  ab  ®  a®  b 

2m  (03  fl-t03/,  )2 


Pa+Pb-2PaPbcos(lab. 


(1) 
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®  b  (P a  Pb  )  <  ®  a®  b  (Pa3~Pb  ) 

03  ab  2m  03  ab  2m 

Therefore,  to  conserve  both  momentum  and  energy  during  the  coalescence  process,  a  quasi-  internal  energy 
is  created,  which  contains  the  randomized  components  of  the  kinetic  energies  of  the  initial  particles.  This  is 
equivalent  to  describing  the  pseudo-particle  with  a  distribution  of  velocities,  such  that  the  temperature  of  an 
equivalent  Maxwellian  distribution  is  :bE  =  Y2kT .  One  could  add  to  the  description  of  the  particle  this 
“thermal”  component,  and  there  are  a  number  of  arguments  in  favor  of  such  a  modification.  A  very 
sophisticated  model  based  on  this  approach  has  been  recently  described  by  Hewett  [11],  and  demonstrated 
high  accuracy  with  fewer  number  of  particles.  Besides  the  reduction  in  computational  cost,  there  is  another 
potential  advantage  in  having  internal  velocity  DF  in  the  pseudo-particle  data  structure.  When  dealing  with 
collisional  plasmas,  the  rates  of  excitation  and  ionization  are  very  sensitive  to  the  energy  distribution 
function  near  the  inelastic  thresholds,  and  can  be  a  critical  part  of  the  overall  plasma  dynamics.  Using  a 
MCC  method  with  only  a  few  particles  may  generate  an  unacceptable  level  of  stochastic  noise;  if  each 
pseudo-particle  also  contains  its  own  velocity  DF,  a  smoother  and  more  accurate  evaluation  of  these  rates 
could  possibly  be  obtained.  Assuming  a  Gaussian  distribution  centered  on  the  mean  particle  velocity  and 
using  a  Taylor  expansion  of  the  cross-section  around  that  velocity,  one  can  easily  evaluate  the  inelastic  rate 
as  function  of  the  pseudo-particle  temperature1)  However,  a  succession  of  merging  processes  would 
eventually  lead  to  excessive  internal  DF  temperatures,  and  assigning  dynamics  based  on  the  mean  velocity 
only  (i.e.  particle  trajectories)  leads  to  inaccuracies.  In  the  absence  of  the  internal  DF  of  each  pseudo¬ 
particle,  it  is  clear  from  the  discussion  above  that  one  cannot  combine  two  particles  into  a  single  one  while 
conserving  mass,  momentum  and  energy.  There  have  been  various  schemes  [19-20]  proposed  to  minimize 
the  energy  error,  but  the  one  described  below  provides  an  exact  conservation  scheme. 


Let  us  first  consider  the  more  general  process  of  combining  an  arbitrarily  large  number  of  particles,  and  to 
consider  the  energy  in  the  center  of  mass  (CM)  frame.  Combining  N  identical  particles,  we  get  a  total  mass 
M  =  y^  Q3;m  ,  and  a  total  momentum  ft  =  p,  .  The  velocity  of  the  center  of  mass  is:  $ CM  =$/  M  . 


Denoting  by  a  prime  the  variables  in  that  CM  reference  frame,  the  individual  constituents  have  a 
momentum :p\  =  pt  -  mU CM  ,  and  contribute  a  kinetic  energy  )2  /  2m  to  the  “thermal”  energy. 


The  latter  can  be  expressed  as  the  difference  between  the  total  energy  contributions  of  all  constituents  and 
the  kinetic  energy  of  the  center  of  mass,  i.e.: 


P2  Hi  P.2 

5 e=  V  to,.  =  y  03 .  Li. 

^  2m  2  M  ^  2m 


(3) 


This  particle  distribution  can  now  be  split  into  two  pseudo-particles  with  half  the  statistical  weight  of  the 
entire  system.  Furthermore,  this  can  be  done  while  transforming  the  thermal  energy  into  the  relative  kinetic 
energy  of  the  two  “children”.  Again,  consider  the  split  in  the  CM  frame,  into  two  pseudo-particles  of  mass 
Ml  2  with  equal  and  opposite  momentum  bp'  (see  Figure  4).  The  total  kinetic  energy  is  now: 


2x- 


(§f')2 


■  =  bE 


(4) 


2(M/2) 

Therefore,  one  can  generate  two  pseudo-particles  with  half  the  mass  and  statistical  weight  of  the  center-of- 
mass,  and  with  equal  and  opposite  momentum  in  the  CM  frame,  with  the  magnitude  is  given  by  (4),  i.e.: 

isf'h  JmEeTi  (5) 


and  with  an  orientation  chosen  at  random  along  the  unit  solid  angle.  Note  that  there  is  no  need  for  a  formal 
transformation  to  the  CM  rest  frame;  the  quantities  of  interest  can  be  evaluated  directly,  using  (3)  for  the 

thermal  energy,  and  with  post-fragmentation  momentums:  jfi±bp' ,  using  (5). 


The  procedure  outlined  produces  two  particles  from  N  particles  without  creating  internal  energies,  and 
while  exactly  conserving  mass,  momentum  and  energy.  As  long  as  N>2,  there  is  an  effective  reduction  in 
the  number  of  particles,  and  while  there  is  no  conservative  process  to  coalesce  2  particles  into  1,  the 
process  of  merging  4  into  2  is  equivalent.  The  reduction  factor  achieved  here  can  take  all  values  (n+2)/2 
starting  from  n=l,  and  is  quite  efficient  for  high  values  of  n. 


1  The  procedures  for  evaluating  the  inelastic  processes  are  described  in  section  3.2. 
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We  conducted  preliminary  tests  of  this  procedure  on  a  simple  planar  diode  case.  Without  particle  merging, 
the  ionization  of  the  gas  rapidly  leads  to  exponential  growth  in  the  number  of  simulated  electrons.  With  an 
aggressive  coalescence  process  in  place,  we  were  able  to  reduce  the  number  of  particles  from  250,000  to 
35,000  with  no  significant  losses  in  accuracy  of  the  results.  Figures  5  and  6  show  the  potential  profiles  and 
particle  charge  densities  obtained  in  both  cases  at  the  same  simulated  time.  Note  that  initially  linear,  the 
potential  profile  has  developed  a  virtual  anode  that  partially  traps  electrons  and  slowly  migrates  towards  the 
cathode  (the  ions  are  immobile  in  this  simulation:  the  apparent  motion  is  a  shift  in  the  ionization  front). 
This  also  causes  the  field  at  the  cathode  to  increase,  and  the  electron  emission  rate  to  increase  as  well.  This 
profile  is  in  agreement  with  analytical  models  [3].  Because  the  potential  is  obtained  as  a  result  of  solving 
Poisson’s  equation,  which  filters-out  short- wavelength  fluctuations,  a  better  test  is  the  comparison  of  the 
charge  densities  themselves.  Although  Figure  8  shows  that  the  statistical  noise  is  slightly  higher,  as 
expected,  in  the  merging  case,  the  results  are  in  excellent  agreement.  We  should  point  out  that  currently 
there  is  no  preferred  selection  of  which  particles  to  be  merged,  and  this  can  results  in  artificial 
thermalization  in  some  cases  (e.g.  beam-beam  interactions);  this  feature  is  currently  being  developed. 


3.2  Collisional-Radiative  MCC 

The  intent  to  implement  a  Collisional-Radiative  model  with  a  PIC  model  for  the  high-energy  plasma 
component  can  lead  to  some  difficulties.  Usually,  collisions  are  treated  within  the  framework  of  the  PIC 
approach  with  a  stochastic  method,  i.e.  Monte-Carlo  Collisions  (MCC).  In  this  approach,  the  probability 
Pc  of  the  modeled  pseudo-particle  having  a  collision  leading  to  some  event,  whether  energy  exchange, 
ionization,  excitation  or  other,  is  compared  to  a  random  number  R  .  If  the  event  is  “probable”  (R  <  Pc  ),  the 
event  is  being  computed,  i.e.  in  the  case  of  ionization,  the  pseudo-particle  looses  the  energy  required  to 
ionize  the  target,  an  additional  electron  is  produced  with  the  same  statistical  weight  as  the  impacting 
particle,  and  the  number  density  of  the  target  is  decreased  by  a  corresponding  amount.  In  the  case  of  a  C-R 
model,  one  must  also  consider  a  multiplicity  of  electronic  levels.  Consider  ionization  from  any  atomic  level 
n  by  electron  impact,  as  shown  in  Figure  7.  In  is  the  ionization  potential  of  that  level,  and  let  us  assume  that 
this  bound-free  transition  leaves  an  ejected  electron  with  a  negligible  kinetic  energy.  The  rate  of  ionization 
due  to  such  transitions  from  level  n  is: 


-^L=Nn\  dzc(z)v-m 

Ot 


(6) 


where  v  =  \l2zlm  (or  v=Jc2  (l+mc2  /£ )”  in  the  relativistic  case),  8  is  the  kinetic  energy  of  the 
electron,  /(e)  is  the  normalized  electron  energy  distribution  function  (EEDF).  After  discretization  of  the 
distribution  function  (see  Figure  8): 


which  we  can  re-write  as: 


=  Nn L ° 


dt 

dNt_ 

dt 


!m 


(7) 


=  N, 


'I*' 


(k)  '  ne( k) 


where  ne(k)  =  Aekf(£k)  =  J <*•/(£) 


(8a) 


and  ki(k)  =a(£]c)- ^e^/m  (8b) 

The  rate  of  change  is  obtained  as  the  sum  of  elementary  processes,  such  that  electrons  from  kth-bin  are 
moved  to  another  energy  bin  m  such  that: 

&k  ~  £/n  =  1  n  (9) 


For  each  elementary  process,  we  have: 


dNi  _ 

d ne(k ) 

dt 


~~dt L  =  +kWneWNn 

=  ~kKk)‘  ne(k)'  Nn 


(10a) 

(10b) 
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dn£ 


e(m) 


dt 


dn, 


e(0 ) 


■-+ki(kyne(ky  Nn 

(10c) 

=  +ki(ky  ™e(ky  Nn 

(lOd) 

7  +  y\fc8r  A 

"  5  Ni  " 

rhs[ 

K  7+jy*  8* 

5  N„ 

rhs„ 

M  O 

k) 

rhse(k) 

1 - 

£ 

Oo 

M 

M 

dt 

The  last  equation  describes  the  production  of  electrons  with  zero  kinetic  energy  (the  ejected  atomic 
electron).  The  interesting  aspect  of  eqs  (10)  is  that  one  can  construct  a  complete  system  of  kinetic  equations 
and  solve  implicitly,  i.e.  when  expanding  the  RHS  to  first  order  in  Sn,§N  : 


(11) 


One  can  then  use  time  scales  larger  than  the  characteristic  time  scale  of  ionization/recombination,  and  by 
extension,  excitation/de-excitation  for  a  full  C-R  model.  However,  in  the  PIC  approach  the  EEDF  is 
approximated  by  a  summation  of  Dirac  delta-functions,  one  for  each  pseudo-particle,  i.e.: 

ne=jdef(e)=jd^g(S)=^Y,mP83-^P)-8(SLp)  (12) 

0  0  P 

where  05  p  is  the  statistical  weight  of  the  particle.  Consider  now  the  standard  MCC  approach.  Suppose  that 
one  is  considering  the  ionization  from  an  excited  level  which  has  a  low  population  density  ( N  n  «  N() ). 
The  rate  of  ionization  due  to  a  given  pseudo-particle  is  given  by: 
dNi 


Oiy  ®  n  / -  05  n 

-L=Nn^Lc(zp)-j2^  =  k;(p)Nn^L 


(13) 


The  probability  of  this  event  to  occur  within  a  time  step  At  is 

Pc  =l.-e-v‘At  ~v cAt  =kKp)NnAt  (14) 

and  is  small  when  the  level  population  is  low.  Nevertheless,  there  will  be  occasional  events  when  ionization 
from  that  level  does  occur.  The  MCC  procedure  would  then  add  to  the  ion  density  and  remove  from  the 
density  of  excited  levels.  However,  the  number  density  of  that  excited  state  would  change  at  that  time  by: 

03  , 

dN„  =- — —  (15) 


P 

V 


which  can  easily  lead  to  negative  population  densities  when  03  p  IV  »  Nn  . 


One  could  imagine  a  modified  procedure  where  conservation  properties  (mass,  energy)  are  conserved  in 
average,  by  implementing  a  small  change  in  the  population  density  at  each  time  step,  irrespective  of 
whether  the  collision  event  is  accepted  or  not.  In  that  case,  the  elementary  change  would  be: 

dNn=AL-kmNnto)  (16) 

Therefore,  if  the  term  in  parentheses  is  sufficiently  small,  the  density  of  excited  states  can  remain  physical. 
However,  this  is  not  always  true;  even  if  Nn  decreases  as  the  level  energy  increases,  the  corresponding 
cross-section  of  ionization  also  increases.  Therefore,  the  product  is  not  always  very  small,  and  there  is  still 
a  chance  that  the  gradual  change  technique,  besides  the  problem  with  conservation  properties,  would  not  be 
able  to  prevent  negative  population  densities.  What  is  needed  is  an  implicit  C-R  method  similar  to  (11), 
which  also  works  when  some  of  the  rates  are  provided  by  interaction  with  pseudo-particles  from  a  PIC 
model.  There  are  several  approaches  to  this  problem,  which  are  currently  being  investigated. 


3.2  Multi-Temperature  Model 

The  simplest  models  for  the  convection  of  the  plasma  fluid  variables  are  finite-volume  based;  in  this 
approach,  the  fluxes  are  evaluated  at  the  cell  interfaces  and  Gauss  law  is  used  to  compute  the  changes  in  the 
cell-averaged  variables.  This  so-called  conservative  form  can  be  written  as: 
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—  +  Ff  =  #  (17) 

dt 

with  Q  the  array  of  conserved  variables,  F  the  array  of  convective  fluxes  and  £$4he  array  of  local  source 
and  exchange  terms.  This  approach  is  well-suited  for  the  arbitrary  grid  topology  described  in  section  <>; 
although  structured  solvers  are  currently  being  implemented,  we  envision  at  a  later  stage  to  supplement  the 
methods  used  to  compute  the  fluxes  by  others  designed  for  unstructured  grids.  Presently,  the  fluid  models 
being  implemented  are  single-  and  multi-temperature  versions  of  a  second-order  Total  Variation 
Diminishing  (TVD)  scheme.  For  convection  processes  only  (i.e.  the  Euler  equations)  and  ignoring  the 
source  terms  due  to  collisional  exchanges  of  mass  (chemistry)  and  energy  (relaxation),  we  have: 


Q  = 


"  pj  ’ 

P_£^a 

0 

p  ua 

0 

M 

= 

M 

M 

E 

u^(E+P) 

Je  _ 

U^Ee 

P  P 
-  pew  K 

(18) 


The  total  energy  E  of  the  plasma  includes  a  purely  thermal  component  E  and  the  kinetic  energy. 
E  =  E+j/2  p  u2  .  The  pressure  tensor  can  be  defined  as  the  sum  of  two  contributions: 

(19) 

static  dynamic 

For  a  two-temperature  plasma  [21],  there  are  two  contributions  to  the  pressure,  and  to  the  plasma  internal 
energy,  i.e.  P  =Ph  +  Fe ,  with  (Ye  =  5/3): 

Ph  -  (ifh  -i)Eh ;  Pe  =  (ye-l)Ee  (20) 

This  implies  that  a  jump  in  pressure  across  a  cell  interface  can  be  expressed  by: 

AP  =<yh  -l)(AE-AEK-Me)  +  (ye  -l)AEe  (21) 

P 2  P 2 

where:  A EK  =  $ •  A(pjP)  -  - l—  Ap  =  pf  -  A$  +  - l—  Ap  (22) 

The  linearization  AF  =  AoAQ  forms  the  basis  of  the  numerical  schemes  based  on  approximate  Riemann 
solvers  [22]  (so-called  Godunov  schemes  [23]),  such  as  the  TVD  scheme  [24].  Notice  that  equation  (18) 
contains  a  source  term  on  the  Right  Hand  Side  (RHS)  for  the  electron  energy  equation.  This  term  is  a 
consequence  of  the  one-fluid  approximation,  and  is  responsible  for  the  adiabaticity  of  the  electron  gas 
component.  Notably,  since  the  electron  gas  is  always  subsonic,  there  is  no  shock  heating  of  the  electrons 
across  a  viscous  shock,  where  there  can  be  a  considerable  jump  in  ion  temperature  [25].  The  complete 
equation  for  the  electron  energy,  combined  with  the  conservation  of  mass  for  the  electrons,  leads  to: 

T?  {D,Te)=(2l  {Dtne) ,  (23) 

where  Dt  =  d  t  +  u  T  is  the  covariant  derivative.  Therefore  the  electron  heating  (cooling)  obtained  as  a 
result  of  compression  (expansion)  is  always  adiabatic :  log(7^)  =  (y  e  -l)log(ne) .  The  single-fluid 

formulation  has  an  important  problem:  the  PeVu  term  on  the  RHS  is  not  an  exact  flux  differential  and 
therefore  cannot  handle  discontinuities.  A  modified  formulation  using  an  entropy-like  variable  [26]  leads  to 
a  simpler  and  more  accurate  treatment.  We  defined 

Pe  .  -  Pe 


=- 


Ae-l 


(24) 


The  convection  of  this  quantity  yields: 


pe  was  used  in  the  original  formulation  [1],  in  which  case  se  was  the  tme  electron  entropy,  but  the  advantage  here  is 
that  the  variables  in  (1)  remain  definite  when  the  ionization  fraction  goes  to  zero  ( Pe  -^o ,  but  p  *0). 
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3* 


=  (ye-I)p1  Y'[a 


tEe  + 


V(f?Ee  )]+(!-' 


Ve) 


(25) 


Using  the  equation  for  the  electron  thermal  energy  dtEe  +  V  •  (uEe  )  =  -pe  (v  •  u),  we  can  easily  verify  that 
the  modified  electron  entropy  satisfies  a  simple  conservation  equation,  i.e.: 

dtSe+V(uSe)  =  0  (26) 

The  TVD  scheme,  modified  to  include  this  conservation  equation,  is  completely  described  for  the  adiabatic 
speed  of  sound  in  the  two-temperature  (ion-acoustic)  plasma: 


a 


2 


(27) 


The  use  of  the  entropy  variable  may  be  convenient  for  the  convective  processes,  but  not  for  electron  heat 
conduction  and  inelastic  energy  exchanges.  The  latter  make  it  more  appropriate  to  retain  the  electron 
thermal  energy  as  the  conserved  variable.  However,  the  operator- splitting  approach  allows  us  to  compute 
the  convective  processes  differently.  Thus,  the  TVD  scheme  can  be  used  as  described  above,  using  the 
conservation  equation  for  Se ,  and  yielding  a  change  in  the  modified  entropy  variable: 
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surfaces 


From  which  a  corresponding  change  in  electron  thermal  energy  due  to  convection  is  computed: 
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ye-i 


-l  s 
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where  5p ,  5 Se  are  the  results  of  convection  only,  and  are  computed  by  the  TVD  scheme.  This  overall 
approach  is  more  accurate  and  uses  less  operations  than  the  previous  scheme  based  on  both  convective  and 
advective  terms;  the  relative  errors  for  the  same  cases  of  shock  and  rarefaction  waves  are  now  of  the  order 
of  10  3  for  a  shock  and  3xl0'5  across  a  rarefaction.  Figure  9  shows  an  example  of  a  test  of  the  multi-fluid 
model,  including  relaxation  effects  and  precursor  heat  conduction. 


3.3  Multi-Fluid  Model 

As  an  alternative  to  a  PIC  model  for  the  energetic  electrons,  one  can  envision  separating  the  electrons  from 
the  rest  of  the  plasma  while  keeping  a  continuum  model.  This  is  the  basis  of  the  multi-fluid  model;  it 
assumes  that  the  electrons  are  sufficiently  collisional  within  themselves  to  maintain  a  DF  that  is  close  to  a 
Maxwellian,  drifting  at  a  mean  velocity  that  is  governed  by  both  inertial  and  electrostatic  effects.  The  main 
advantage  is  that  computations  can  be  much  faster  than  with  a  PIC  model,  while  still  being  able  to  account 
for  space-charge  effects.  The  disadvantage  is  that  the  time  scale  is  governed  by  a  combination  of  the  mean 
electron  velocity  and  the  speed  of  sound  in  the  electron  gas.  Therefore,  for  this  approach  to  be  useful,  an 
implicit  fluid  solver  must  be  devised  for  the  electrons,  coupled  to  the  electric  field.  A  first  attempt  at 
devising  such  a  scheme  has  been  made.  Described  in  more  details  in  Appendix  A,  two  versions  of  the 
scheme  have  been  designed;  one  which  entails  the  inversion  of  7x7  matrices,  and  is  linear-order  in  the  time 
step  bt ;  a  second  version  that  uses  4x4  matrices,  but  is  a  second-order  expansion,  i.e.  includes  terms 

proportional  to  5 12 .  Both  algorithms,  i.e.  the  7x7  implementation  given  by  (A20)  and  the  second-order 
expansion  algorithm  described  in  (A22)  were  implemented  and  tested.  The  test  case  consisted  of  a  simple 
one-dimensional  plasma  with  reflecting  boundaries  (“plasma  in  a  box”);  the  initial  state  was  a  fully  ionized 
plasma  without  recombinations,  and  with  all  the  electrons  moved  to  the  left  half  of  the  box.  This  generated 
a  large  internal  electric  field  which  accelerated  the  electrons  towards  the  right  side.  Without  dissipation,  the 
electron  motion  leads  to  an  oscillatory  behavior,  as  the  electrons  accumulate  on  the  right  side,  create  a 
restoring  field  that  accelerates  back  the  other  way,  and  so  on.  This  test  case  is  similar  to  the  one  studied  in 
[Shumlak],  but  with  a  much  more  severe  initial  condition;  the  initial  state  in  [Shumlak]  was  a  slight 
perturbation  of  the  space  charge,  in  order  to  have  linear  superposition  of  the  natural  (Langmuir)  plasma 
oscillations  with  the  background  fluid  at  rest.  In  our  case,  the  plasma  oscillations  have  a  much  larger 
amplitude,  and  non-linear  effects  (steepening  of  the  profiles)  can  be  detected;  nevertheless,  it  will  be  seen 
that  one  can  still  obtain  an  excellent  agreement  with  theory.  This  test  case  is  also  a  much  more  severe  test 


12 


35th  AIAA  Plasmadynamics  and  Lasers  Conf. 


AIAA  2004-2166 


of  the  stability  of  the  algorithms.  Indeed,  it  was  found  that  the  7x7  algorithm  (A20)  would  collapse  when 
the  time  step  was  too  large  compared  to  the  plasma  frequency  or  the  Courant  condition  (\u\+c)bt  <  Ax .  The 
second-order  expansion  scheme  provided  greater  stability,  but  damped  the  oscillations  when  the  time  step 
was  such  that  (\u\+c)bt  >  Ax ,  co pebt  »1 .  This  is  to  be  expected,  since  for  large  time  steps  the  algorithms 
cease  to  be  time-accurate.  It  is  not  yet  exactly  clear  why  this  damping  occurs  (albeit  at  a  slow  rate)  for  time 
steps  just  below  the  Courant  condition.  Further  investigations  are  planned  for  the  near-future. 


3.4  MHD 


An  ideal  MHD  solver  is  also  in  the  preliminary  stages  of  being  implemented.  Based  on  previous 
development  work,  the  solver  is  shock-capturing  and  has  been  successfully  tested  on  the  standard  Brio  & 
Wu  problem,  both  with  aligned  and  misaligned  (rotated  problem)  grids.  The  scheme  is  based  on  a  TVD 
scheme  with  the  MHD  characteristic  waves,  and  uses  flux-constrained  transport  concepts  [27]  to  maintain 
the  divergence  of  the  B -field  down  to  machine  precision  levels.  The  method  works  well  on  body-fitted 
grids,  and  is  being  extended  to  include  multi-temperature  and  resistive  effects.  One  of  the  key  problems,  for 
example,  is  to  have  the  scheme  valid  in  the  limit  of  either  small  or  large  conductivity.  The  ideal  MHD 
equations  are  then  obtained  by  transforming  the  expressions  of  the  Lorentz  force  and  Joule  heating  into 
modified  fluxes  for  conservation  laws  of  combined  plasma  and  field  quantities.  Resistive  effects  appear  in 
the  dynamics  of  the  field  and  total  energy,  i.e.: 


dt 


[«|3Ba-«aB|3]=V|3 

_L(VPBa  -V01#13 ) 

(30a) 


2  [Ep +Em  ]+  V  p  [  u  P  (Ep  +Em  )+ua(P “P  +T$  )]=  Vc 


dt 
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ap0 
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where  em  =  B 1  /  2fi0  is  the  magnetic  energy,  the  pressure  tensor  Pap  is  defined  by  (19),  and 

t  oc(3  _  B  s  /oi  \ 

1  M  - - °a(3  “ -  (31) 

2\i0  Fo 

is  the  Maxwell  tensor.  When  the  conductivity  is  small,  the  magnetic  diffusion  and  associated  resistive 
heating  (right-hand- side  terms  in  30)  become  very  large.  These  must  be  computed  with  an  implicit  method 
in  order  to  prevent  arbitrarily  small  time-step  limitations.  A  test  of  this  magnetic  diffusion  scheme  is  shown 
in  Figure  10.  This  test  consists  of  a  high  pressure  core  region  where  the  field  is  also  inverted.  The  field 
discontinuity  is  being  diffused  away,  and  as  the  magnetic  field  decays,  its  energy  is  converted  into  plasma 
heating  (here  the  convective  terms  are  not  being  computed).  The  test  compared  a  time-accurate  explicit 
scheme  with  the  implicit  version  at  various  time  steps;  the  results  show  that  despite  the  large  time  steps  the 
implicit  scheme  is  very  accurate,  even  at  the  initial  discontinuity  where  the  rate  of  heating  is  the  largest. 

Comparisons  between  the  loosely-coupled  approach  ( jxB ,  j  •  E )  and  the  full  MHD  solver  (eqs.  30)  were 
also  successfully  performed.  This  MHD  solver  will  be  upgraded  to  include  approaches  for  both  low  and 
high-beta  plasmas,  and  incorporated  into  the  general  hybrid  model. 


4 .  Future  W ork  and  Conclusions 

Other  solvers,  not  described  here,  have  been  or  are  being  implemented  (for  example,  a  GMRES  method  for 
Poisson’s  equation),  and  discussions  of  the  Equation-Of-State  (EOS)  and  thermo-chemical  database  have 
been  set  aside.  More  importantly,  issues  of  multiple  time-scales  have  not  been  described;  these  are  readily 
apparent  when  considering  the  dynamics  of  both  the  bulk  plasma  and  the  energetic  electrons.  The  latter 
evolve  on  much  faster  time  scales  than  the  ions  and  neutrals,  which  can  be  considered  as  “frozen”  during 
the  time  scale  5 1  of  transport  of  the  electrons  in  the  PIC  model.  It  is  not  desirable  to  completely  eliminate 
the  PIC  approach  in  favor,  for  example,  of  a  ballistic  model  -  where  electron  trajectories  are  computed  one 
by  one  in  a  given  field  configuration  -  since  this  model  would  not  accurately  account  for  fast  transient 
effects  (e.g.  beam  instabilities).  However,  it  is  necessary  to  being  able  to  sub-cycle  the  PIC  model,  such 
that  other  components  are  updated  on  a  time-scale  At  that  is  governed  by  their  own  transport,  or  by  the 
coupling  with  the  fast  (PIC)  component.  For  example,  while  the  PIC  model  is  being  sub-cycled,  the  rates  of 
ionization  and  excitation  must  be  integrated  during  that  time,  until  the  projected  changes  in  population 
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levels  become  significant.  This  is  not,  however,  a  trivial  matter;  during  the  ionization  process,  new  free 
electrons  can  themselves  become  accelerated  by  the  field  and  become  part  of  the  PIC  component.  Thus, 
during  time  scales  that  can  be  much  smaller  than  characteristic  times  of  heavy  particle  transport  or  relative 
changes  in  composition,  the  inelastic  coupling  can  affect  the  time-dependent  dynamics  of  the  PIC 
component,  and  must  be  considered.  The  impact  on  overall  scheduling  and  parallelization  strategies  must 
still  be  evaluated.  Therefore,  while  small  time  steps  can  be  used  for  testing  and  validation  in  the  near  future, 
an  effective  strategy  for  dealing  with  the  multiple  time  scales  will  be  the  next  major  development  task. 

In  conclusion,  the  development  of  a  general-purpose  hybrid  model  is  well  under-way,  and  is  rapidly 
building  up  capabilities.  The  strong  Object-Oriented  approach  is  absolutely  necessary  to  the  success  of  this 
project,  and  the  use  of  a  Java-based  top-level  architecture  provides  an  elegant  solution  to  the  problems  of 
portability  and  accessibility  while  maintaining  a  minimum  level  of  performance;  the  approach  also  lends 
itself  to  the  modem  paradigm  of  grid  computing,  allowing  remote  steering  and  flexible  access  to  computing 
services.  Several  numerical  issues  specific  to  the  cases  of  interest,  namely  highly  energetic  particles  inside 
a  dense,  collisional  bulk  plasma,  have  been  discussed.  New  models  and  approaches  to  grid,  particle  and 
fluid  data-stmctures,  particle  transport,  and  kinetics  of  particle-fluid  interaction  are  being  developed.  A 
wide  variety  of  fluid  models  is  also  being  implemented  and  tested.  The  resulting  capability,  along  with 
upcoming  implementations  of  a  relativistic,  fully  electro-magnetic  PIC  and  Maxwell  solvers,  is  projected  to 
be  very  useful  in  the  study  of  the  fundamental  physics  of  new  classes  of  plasmas  discharges. 
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Appendix  A:  Electron-Fluid  Model 


Consider  the  conservation  equations  for  density  and  momentum  only  in  one  -direction: 

dtn+dx(nu)=0  (A-l) 

dt(nu)  +  dx(nu2  +  pe/me)  =0  (A-2) 

define  an  eigenvalue:  c2e  -  kTe  / me  which  is  the  speed  of  sound  in  an  isothermal  gas.  Then  (A-l)  and  (A- 
2)  can  be  written  as  dtQ  +  dxF  =  0  ,  where  Q  is  the  vector  of  conserved  variables  and  F  is  the  flux  in 
conservative  form: 


3, 


We  define  the  transformation  matrices: 


X  = 


+  d. 


n{u2+c l) 


=  0 


and  the  diagonal  eigenvalue  matrix: 
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The  Jacobian  of  the  system  (A-3)  is  obtained  as:  A  =  X  AX.  The  characteristic  jumps  are: 

a  =XAQ  = 


'  tzl  km  A(nM) 

2c  f  2c 
A  (nu) 
2c 


(A-3) 
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In  order  to  construct  a  Godunov  scheme  that  uses  flux-limiters  on  the  characteristic  variables,  it  is 
necessary  to  verify  that:  (a)  A Q  =  . 

For  the  second,  we  first  note  that: 


necessary  to  verify  that:  (a)  AQ  =  X  1  a  ;  and  (b)  A F  =  X  1  •  A  •  a  .  The  first  relation  is  trivial  to  verify. 


Aa  = 


Then,  after  some  algebra: 
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X  1  Aoc  m 


A  (nu) 

-  (u  2-c2)An  +  2uA(nu) 
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A(nu2)+c2An 
2 


(A-7) 


The  second  term  is  equal  to  the  flux  component  of  (A-3)  as  long  as  c  is  invariant,  in  which  case 

c2 An  =  A {nc2 ) .  This  is  true  because  the  system  (A- 1,2)  describes  an  isothermal  gas:  the  pressure  is 
controlled  by  the  temperature,  which  is  an  external  variable  here:  there  is  no  equation  of  state  (EOS)  that 
determines  the  temperature  from  only  n,  u ,  and  there  is  no  conservation  equation  for  the  energy  or 
pressure.  The  system  can  be  generalized  to  more  than  one  dimension,  where  (A-2)  extends  to: 

dt(nua)  +  d  $  ( nuau ^  +  8°^  pe/me)  =  0  (A-8) 

The  system  (A-3)  becomes: 
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(A-9) 


The  flux  along  a  direction  n  (with  t,s  the  correspondingly  rotated  unit  vectors  normal  to  n)  can  be 
obtained  through  the  following  matrix  definitions: 


X  = 


and  the  eigenvalue  matrix  is: 
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With  these  matrices,  the  flux-limited  solver  can  be  constructed.  An  implicit  algorithm  can  also  be 
constructed,  such  that: 
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(A- 11) 


where  A+,A  are  the  Jacobians  for  positive  and  negative  eigenvalues  only.  They  are  constructed  by 
=X~1-A±-X  ,  with:  A+  =  {K  ,max(<9,?ia),K } ,  and  A“  =  {K  ,min(0,Xa),K  } . 

We  can  now  turn  our  attention  to  the  right-hand  sides  of  the  conservation  equations,  which  have  been 
ignored  so  far.  We  will  consider  first  the  electrostatic  approximation,  which  implies  that  only  an  electric 
field  is  present,  and  all  magnetic  fields  can  be  neglected.  The  conservation  equation  (A-l)  does  not  change, 
but  the  momentum  equation  is  now: 
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is  the  frictional  force  and  v  en  ,v  ei  are  the  collisional  frequencies  for  momentum  exchange.  The  latter  can 

be  simplified  when  the  ions  and  neutrals  are  both  at  the  same  velocity  U  ;  if  the  plasma  is  at  sufficiently 
low  density  for  an  ion-slip  to  occur,  the  approach  described  below  can  easily  be  extended  to  the  more 
general  case.  For  simplification,  we  also  neglect  the  thermal  force,  which  can  also  be  easily  added  later  on. 

The  second  term  on  the  left-hand- side  is  the  advection  term  which  was  treated  in  the  previous  section.  The 
contributions  of  each  term  on  the  RHS  can  be  examined  one  by  one.  The  easiest  is  the  collisional 
momentum  exchange.  One  can  write  this  contribution  as  follows: 

8(ne^)=8tveh[neUa -(«,«“)] 

Evaluating  terms  on  the  RHS  at  the  advanced  time  level  will  yield: 

8  (neUe  )  =5 tv eh  \neUa-(neu%  )]  + (htv  ehUa  )hne  -  (htv eh)h(nei%  ) 

The  term  involving  the  electric  field  can  be  written  in  a  similar  form, 


8  (neu®  )  =  -ht-^—neEa  -ht-^—nehEa  -St  Eahne 
me  me  me 

and  we  need  only  to  express  the  change  in  the  electric  field.  Using  one  of  Maxwell’s  equations: 

p 


(A- 14) 
(A- 15) 
(A- 16) 


£„ - +  J  =  - 

°  dt 


and  neglecting  the  magnetic  field  (electrostatic  approximation),  this  leads  to: 

Z0dtEa  =-e[ziniUf  -neu%) 

Again,  expressing  the  RHS  at  the  advanced  time  level  yields: 

8 Ea  =  — —St  {zinjUfJ  - neu “  )+  —hth(neu®) 

£0 

One  can  therefore  construct  a  point-implicit  solution  of  the  form: 


- 

7 

0 

0 

(hn  ^ 

(  \ 

0 

eEPht 

+ 

me 

~^ehUa 

MtVeh 

enpht 
+  e 

me 

h(neu “) 

8  Ea 

v  ) 

= 

8 r  v  ehne  (Ua  )  -8 t^-Ea 

me 

- 

0 

-  eh t/£0 

1 
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(A- 17) 

(A- 18) 
(A- 19) 

(A-20) 


Note  that  this  is  a  point-implicit  part  due  to  the  Taylor  expansion  of  the  RHS.  One  should  also  include  the 
advection  process,  i.e.  the  matrices  defined  in  (A-l  1).  Without  the  energy  conservation  equation,  this  would 
lead  to  a  7x7  matrix  system:  3(momentum)  +  3(electric  field)  +  1  (density).  However,  since  the  evolution  of 
the  electric  field  does  not  involve  neighboring  points,  it  would  be  simpler  to  solve  for  the  change  in  electric 
field,  and  include  this  solution  into  the  evolution  of  the  momentum.  Inserting  (A- 19)  into  (A  16),  we  obtain: 

[7+  (tipe§t2] -h{neu^  )  =  ~(eht/me)neEa  +co 2pe§t2  - neu J  )-(eht/me)Ea 8ne  (A-21) 


and  using  J  =  e(ZiniUi  -neUe  ),(A20)  becomes: 


7  0 

+  ( ebt/me )Ea  -8 tv ehUa  TfV  ^Sf-KD^St2 


hne 

h(neu^) 


-  (eht/me  )neEa +us2eht2  \Ja  I e  J 


(A-22) 
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Cathode  Cathode  Fall  Negative  Glow 

Figure  1 :  Schematic  of  physical  processes  near  cathode. 


Figure  2 :  Possible  configurations  for  the  Client/Master/Server  architecture 


Block  1  Block  2 

Figure  3:  Block  communication  using  source  and  target  patches.  The  dashed 
arrows  represent  referencing  while  the  solid  arrows  signify  remote  method  calls. 
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Figure  4:  Schematic  of  fragmentation  process  in  Center  of  Mass  frame  after  coalescence 
of  many  individual  pseudo-particles. 


Figure  5:  Electric  potential  in  diode  10  ns  after  start,  with  and  without  particle  coalescence. 


Figure  6:  Particle  (electron)  and  fluid  (ions)  charge  densities  10  ns  after  start  of  computations, 
with  and  without  pseudo-particle  coalescence.  Reduced  statistical  fluctuations  in  the  former  case 
indicate  the  approximately  9-times  larger  number  of  simulated  particles. 
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continuum 


n=l 


n-0 


Figure  7:  Schematic  of  ionization  scheme  in  C-R  model. 


Cell  Number 

Figure  9:  Temperature  profiles  behind  strong  shock  with  relaxation  and  electron  heat  conduction. 
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?  pe_1  =  226.4  ?  0.5  ps  (theor.  =  226.7  ps) 


Figure  10:  Plasma  oscillations  (electron  velocity  at  center  of  ID  box)  for  explicit  and  implicit  schemes. 
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Figure  11:  comparison  of  heating  rates  (pressure  rise)  for  explicit  and  implicit  magnetic  diffusion  schemes. 
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